#include "ScreeningMola.h"


void ScreeningMola::iniciaMatriz( REAL max_rad, REAL min_r[2], REAL max_r[2] ) {

        assert(Cell==NULL);

        MinR[0] = min_r[0];
        MinR[1] = min_r[1];
        MaxR[0] = max_r[0];
        MaxR[1] = max_r[1];
        CellSize = max_rad;
        CellSize *= 8.0;
        Rows = 1 + ((MaxR[1] - MinR[1]) / CellSize);
        Cols = 1 + ((MaxR[0] - MinR[0]) / CellSize);
        printf ("\tScreening -> Matriz %dx%d\n\t",Rows,Cols);
        int qtdEsferas=0;
        assert (Spheres!=NULL);
        ImprimeMem("Screening",((Rows+2)*(Cols+2)*sizeof(SPHERE*)));
        Cell = new SPHERE**[Rows+2];
        for (int i = 0; i < Rows+2; i++) Cell[i] = new SPHERE*[Cols+2];
        #ifdef MOLAS_SEMAFORO
        for (SPRING* s = Springs->Begin; s != NULL;s = s->Next) omp_init_lock(&(s->semaforo));
        molasDistantes = new vector<SPRING*>[omp_get_num_procs()];
        #endif
}
void ScreeningMola::CriaMatriz()
{
    assert(  Spheres != NULL );

    for (int i = 0; i < Rows+2; i++)
    {
        for (int j = 0; j < Cols+2; j++)
        {
            Cell[i][j]=NULL;
        }
    }
    //Cria a matriz esparsa do screening
    for (SPHERE* s = Spheres->Begin; s != NULL; s = s->Next) {
        long row = Max( Min( 1 + ( (s->R[1] - MinR[1]) / CellSize ), Rows-1 ), 1 );
        long col = Max( Min( 1 + ( (s->R[0] - MinR[0]) / CellSize ), Cols-1 ), 1 );

        s->CellNext = Cell[row][col];
        Cell[row][col] = s;
    }
}

void ScreeningMola::CriaMatrizParalelo()
{
    #pragma omp for
    for (int i = 0; i < Rows+2; i++)
    {
        for (int j = 0; j < Cols+2; j++)
        {
            Cell[i][j]=NULL;
        }
    }
    //barreira implícita
    #pragma omp single
    {
        //Cria a matriz esparsa do screening
        for (SPHERE* s = Spheres->Begin; s != NULL; s = s->Next) {
            long row = Max( Min( 1 + ( (s->R[1] - MinR[1]) / CellSize ), Rows-1 ), 1 );
            long col = Max( Min( 1 + ( (s->R[0] - MinR[0]) / CellSize ), Cols-1 ), 1 );

            s->CellNext = Cell[row][col];
            Cell[row][col] = s;
        }
    }
    #pragma omp barrier
}


void ScreeningMola::calcular()
{
    #pragma omp single
    {
        if (Springs->Begin!=NULL) valorDaVez=!valorDaVez;
        //CriaMatriz();
    }
    CriaMatrizParalelo();
    if (Rows>Cols)
    {
        #pragma omp for schedule (dynamic,1)
        for (int i = 1; i <= Rows; i++)
        {
            for (int j = 1; j <= Cols; j++)
            {
                //Pego cada elemento de uma célula
                for (SPHERE* s = Cell[i][j]; s != NULL; s = s->CellNext)
                {
                    checkSprings(s,i,j);
                }
            }
        }
    }
    else
    {
        #pragma omp for schedule (dynamic,1)
        for (int j = 1; j <= Cols; j++)
        {
            for (int i = 1; i <= Rows; i++)
            {
                //Pego cada elemento de uma célula
                for (SPHERE* s = Cell[i][j]; s != NULL; s = s->CellNext)
                {
                    checkSprings(s,i,j);
                }
            }
        }
    }
    #pragma omp single
    {
        ProcessaDistantes();
    }

}

void ScreeningMola::ProcessaDistantes()
{

    #ifndef MOLAS_SEMAFORO
    for (int x=0; x != molasDistantes.size(); x++ )
    #endif

    #ifdef MOLAS_SEMAFORO

    for (int x=0; x != omp_get_num_procs(); x++ )

        for (int y=0; y != molasDistantes[x].size(); y++ )
        #endif
        {

            #ifndef MOLAS_SEMAFORO
            SPRING *k = molasDistantes[x];
            #endif

            #ifdef MOLAS_SEMAFORO
            SPRING *k = molasDistantes[x][y];
            #endif
            REAL Dist, Dist2, F;
            REAL R[2], V[2];

            k->Translate();

            R[0] = k->R2[0] - k->R1[0];
            R[1] = k->R2[1] - k->R1[1];
            Dist2 = DOT( R, R );
            Dist = sqrt(Dist2);

            V[0] = k->S2->V[0] - k->S1->V[0];
            V[1] = k->S2->V[1] - k->S1->V[1];
            F = ( k->K * (k->L0 * Dist - Dist2) - k->C * DOT( V, R ) ) / Dist2;

            k->S1->Force( -F, R );
            k->S2->Force( F, R );

            k->S1->Torque( -F * (PERP_DOT( R, k->R1 ) - PERP_DOT( R, k->S1->R )) );
            k->S2->Torque( F * (PERP_DOT( R, k->R2 ) - PERP_DOT( R, k->S2->R )) );

            k->MeasureStress( Dist );
            if (k->Stress > k->Mat->Sigma_E) k->L0 += k->L0 * (k->Stress - k->Mat->Sigma_E) / k->Mat->Young;
            else if (-k->Stress > k->Mat->Sigma_E) k->L0 += k->L0 * (k->Stress + k->Mat->Sigma_E) / k->Mat->Young;

            if (Dist > k->FractureT || Dist < k->FractureC)
            {
                md->Delete(k);
            }
        }

    #ifndef MOLAS_SEMAFORO
    molasDistantes.clear();
    #endif
    #ifdef MOLAS_SEMAFORO
    for (int x=0; x != omp_get_num_procs(); x++ )
            molasDistantes[x].clear();
    #endif

}
void ScreeningMola::checkSprings(SPHERE *s,int l, int c)
{
    for (SPRING* k = s->Springs; k != NULL; )
    {
        SPRING* prox=k->SNext(s);
        int row,col;
        if (s!=k->S1)
        {
            row = Max( Min( 1 + ( (k->S1->R[1] - MinR[1]) / CellSize ), Rows-1 ), 1 );
            col = Max( Min( 1 + ( (k->S1->R[0] - MinR[0]) / CellSize ), Cols-1 ), 1 );
        }
        else
        {
            row = Max( Min( 1 + ( (k->S2->R[1] - MinR[1]) / CellSize ), Rows-1 ), 1 );
            col = Max( Min( 1 + ( (k->S2->R[0] - MinR[0]) / CellSize ), Cols-1 ), 1 );
        }
        if ((row!=l)||(col!=c))
        {
                #ifndef MOLAS_SEMAFORO
                #pragma omp critical
                {
                    if((k->foiAcessado)==valorDaVez)
                    {
                        //Evitar que outro processador coloque a mola duplicada na lista
                        k->foiAcessado=!valorDaVez;
                        molasDistantes.push_back(k);
                    }
                }
                #endif

                #ifdef MOLAS_SEMAFORO
                omp_set_lock(&(k->semaforo));
                if((k->foiAcessado)==valorDaVez)
                {
                    //Evitar que outro processador coloque a mola duplicada na lista
                    k->foiAcessado=!valorDaVez;
                    molasDistantes[omp_get_thread_num()].push_back(k);
                }
                omp_unset_lock(&(k->semaforo));
                #endif
        }

        //Evitar computar a mesma mola quando encontrar a outra esfera conectada
        if ((k->foiAcessado)==valorDaVez)
        {
            //printf("real ");
            REAL Dist, Dist2, F;
            REAL R[2], V[2];

            k->foiAcessado=!valorDaVez;
            k->Translate();

            R[0] = k->R2[0] - k->R1[0];
            R[1] = k->R2[1] - k->R1[1];
            Dist2 = DOT( R, R );
            Dist = sqrt(Dist2);

            V[0] = k->S2->V[0] - k->S1->V[0];
            V[1] = k->S2->V[1] - k->S1->V[1];
            F = ( k->K * (k->L0 * Dist - Dist2) - k->C * DOT( V, R ) ) / Dist2;

            k->S1->Force( -F, R );
            k->S2->Force( F, R );

            k->S1->Torque( -F * (PERP_DOT( R, k->R1 ) - PERP_DOT( R, k->S1->R )) );
            k->S2->Torque( F * (PERP_DOT( R, k->R2 ) - PERP_DOT( R, k->S2->R )) );

            k->MeasureStress( Dist );
            if (k->Stress > k->Mat->Sigma_E) k->L0 += k->L0 * (k->Stress - k->Mat->Sigma_E) / k->Mat->Young;
            else if (-k->Stress > k->Mat->Sigma_E) k->L0 += k->L0 * (k->Stress + k->Mat->Sigma_E) / k->Mat->Young;

            if (Dist > k->FractureT || Dist < k->FractureC)
            {
                k=md->Delete(k);
            }
        }
        k=prox;
    }
}
